Merging of globular clusters within inner galactic regions. 
I. Do they survive the tidal interaction? 
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ABSTRACT 

The main topic of this paper is the investigation of the modes of interaction of globular 



clusters (GCs) moving in the inner part of a galaxy. This is tackled by means of high-resolution 
iV-body simulations, whose first results are presented in this article. Our simulations dealt with 
primordial very massive (order of 10 7 M Q ) GCs that were able to decay, because of dynamical 
friction, into the inner regions of triaxial galaxies on a time much shorter than their internal 
relaxation time. To check the disruptive roles of both tidal forces and GC-GC collisions, their 
effects were maximised by considering clusters on quasi-radial orbits and choosing the initial 
conditions so as to give head-on collisions at each passage through the center. 

The available CPU resources allowed us to simulate clusters with different structural param- 
eters and to follow them on quasi-radial orbits during 8 passages across the center. The main 
findings are: i) clusters with an initial high enough King concentration parameter (c > 1.2), pre- 
serve up to 50% of their initial mass; ii) the inner density distribution of the survived clusters keep 
a King model profile; iii) GC-GC collisions have a negligible effect with respect to that caused 
by the passage through the galactic center; iv) the orbital energy dissipation due to the tidal 
interaction is of the same order of that caused by dynamical friction; v) complex sub-structures 
like "ripples" and "clumps" formed, as observed around real clusters. These findings support 
the validity of the hypothesis of merging of GCs in the galactic central region, with modes that 
deserve further careful investigations. 

Subject headings: globular clusters: general, galaxies: kinematics and dynamics, methods: N-body 
simulations 



1 



1. Introduction 

The study of the detailed structure and evo- 
lution of globular clusters (GCs) in galaxies is a 
modern astrophysical concern, which is particu- 
larly suitable for HST (Kundu et al. 1999) and 
for large ground-based telescopes (Gomez et al. 
2001). Some interesting general considerations on 
the structural properties of GCSs in galaxies are 
found in Ashman & Zepf (1997), where also a sum- 
mary of simplified models and results by other au- 
thors on the topic of formation and destruction 
of globular clusters in galaxies is present. While 
the amount of data for GCSs in galaxies is now 
large and rapidly increasing (see, e.g. Davidge & 
van den Bergh 2005; Dirsch et al. 2005; Forbes et 
al. 2004; Fort et al. 1986; Gomez & Richtler 2004; 
Kundu et al. 2004; Olscn et al. 2004; Zepf et al. 
2004), from a theoretical point of view the study 
of the evolution of GCSs in galaxies has been tack- 
led in different ways but it still lacks of definitive 
conclusions. We just recall here the studies by 
Aguilar et al. (1988), Gncdin & Ostriker (1997) 
and Vesperini (2001), which aim at achieving a 
deeper understanding of the effects of the interac- 
tion of clusters with the galactic field, and the in- 
teresting results obtained by Gnedin et al. (1999) 
on the survival of clusters in dependence on their 
initial conditions. 

It is known that, in elliptical galaxies, the most 
important processes causing the evolution of GCs 
are i) the tidal interaction with the global field 
(Murali & Weinberg 1997a), that has also ii) the 
effect of accelerating the relaxation evaporation 
(Murali & Weinberg 1997b); iii) the impulsive in- 
teraction with a compact object in the galactic 
center (on much smaller length-scales), see Os- 
triker et al. (1989), Capuzzo Dolcetta (1993), Ca- 
puzzo Dolcetta & Tesseri (1997), Capuzzo Dol- 
cetta & Tesseri (1999); iv) the dynamical fric- 
tion (df) due to the bulge-halo matter as firstly 
shown by Tremaine et al. (1975) and then gener- 
alised to the triaxial galaxies case by Pesce et al. 
(1992), Capuzzo Dolcetta (1993), Capuzzo Dol- 
cetta & Vicari (2005), who ascertained that df 
plays a much more important role than what pre- 
viously believed; see also Colpi et al. (1999), van 
den Bosch et al. (1999), Penarrubia et al. (2002). 

Actually, Pesce et al. (1992) (hereafter PCV), 
Capuzzo Dolcetta (1993) and Capuzzo Dolcetta 



& Vicari (2005) showed quantitatively (by mean 
of numerical integration of a large set of orbits of 
GCs in self-consistent galactic models) that mas- 
sive enough clusters have short orbital decaying 
times. More specifically Eqs. (Al) and (A2) in 
Capuzzo Dolcetta (1993) indicate the (almost) to- 
tal loss of the orbital energy in less than 5 Gyr 
for 10 6 M Q GCs on low angular momentum or- 
bits having an initial orbital binding energy per 
unit mass E = 0.6, which is (in units of the 
galactic total potential depth) the energy corre- 
sponding to the background stellar velocity disper- 
sion in the assumed Schwarzschild's triaxial model 
(Schwarzschild 1979); this decay time, inversely 
scaling with GC mass, is less than 500 Myr for 
globulars more massive than 10 7 M Q . Inciden- 
tally, the apocentric distance of a quasi-radial or- 
bit (thin box) of E = 0.6 is 6.3r;, while the radius 
of the quasi-circular orbit of same energy is 3.8r{,, 
where is the core radius of the Schwarzschild's 
(1979) model spherical component; the assump- 
tion (as in Schwarzschild 1979) n, = 200 pc gives 
1.3 kpc and 760 pc for the apocenter and circu- 
lar radius, respectively. The more recent paper 
Capuzzo Dolcetta & Vicari (2005) represents a 
generalization and widening of the previous PCV 
and Capuzzo Dolcetta (1993) works and refers to 
a set of model for triaxial galaxies with a density 
core and various axial ratios. For the sake of com- 
parison with the above cited results based on the 
Schwarzschild's model, a semimajor axis of 300 pc 
and the same axial ratios 2:1.25:1 of that model are 
used in the df decay time fitting formulas given in 
Capuzzo Dolcetta & Vicari (2005) to get, for the 
same 10 6 M Q GC of the same initial orbital energy 
and angular momentum, a decay time 10 9 yr, i.e. 
a factor 5 shorter than the one evaluated as de- 
scribed above. Being the df time easily evaluated 
for any GC mass, due to its inverse proportional- 
ity to the mass, we get a decay time ranging from 
10 8 yr to 20 Gyr for an interval of GC masses ex- 
tending from 5 x 10 4 M Q to 10 7 M , implying a 
significant evolution of an initial GC population 
by df, as actually shown by Capuzzo Dolcetta & 
Vicari (2005). 

One application of the present work is to the 
study of the validity of an intriguing scenario 
elsewhere presented (Capuzzo Dolcetta 1993) and 
that naturally emerges from all these considera- 
tions: the high efficiency of the df implies that 
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many massive (M > 5 x 10 6 M Q ) globular clusters 
decay rapidly to the inner galactic region where 
they interact closely with each other and with 
the galactic nucleus, eventually forming, through 
a merging process, a dense and massive super- 
cluster. This scenario, that requires an accurate 
modeling, has important implications both on the 
accretion of massive objects in the galactic cen- 
ter and on their radiative emission as AGN (see 
the discussion in Capuzzo Dolcetta 1993; Capuzzo 
Dolcetta & Vicari 2005). 

One crucial point for this scheme to work is the 
existence of few primordial massive GCs. At this 
regard, we remind that the real problem was, once, 
to understand why very massive GCs were not ob- 
served, indeed. This because the typical Jeans 
mass in a primordial, virialized, gaseous galac- 
tic halo (T ~ 10 5 + 10 6 K) is > 10 7 M , which 
is too large respect to the observed (Milky Way) 
GC mass values. This led to the theoretical at- 
tempt to understand why such massive GCs were 
not observed, on the basis of various schemes and 
scenarios of GC formation (Fall & Rees (1985), Vi- 
etri & Pesce (1995)) mainly invoking the emcience 
of cooling mechanisms reducing the temperature 
to low enough values. The new, relevant, point 
is that there are now quite a few recent papers 
showing, indeed, the existence of young massive 
clusters in Antennae (Fritze-v. Alvensleben 1999), 
in the Magellanic Clouds, M33, Fornax dSph (dc 
Grijs et al. 2005). Young massive clusters have 
been found, also, in M31 (Fusi Pecci et al. 2005); in 
this galaxy, Huxor et al. (2005) discovered bright 
clusters with anomalously large half-mass radii 
(w 30 pc). 

In particular, Fritze-v. Alvensleben (1999) 
showed that the mass function of young clusters 
in the Antennae extends up to few times 10 7 Mq . 
The existence of very massive GCs does not seem 
to be a peculiarity of young systems, only. Ac- 
tually, Harris & Pudritz (1994) give the evidence 
of the presence of 10 7 M in the giant ellipti- 
cal M87 as well as in Virgo ellipticals, being the 
cumulative (sum over 3 ellipticals in the cluster) 
mass function of GCs extended up to high masses. 
The observed presence of very massive GCs fits 
with the recent theoretical-numerical findings by 
Kravtsov & Gncdin (2005) who deduce a depen- 
dence of the most massive GC mass (M max ) on 
the parent galaxy mass (M g ), M max oc M^ 29 , 



such that M max ~ 10 7 M for M g ~ 2.6 x 10 11 
M . Moreover, they also find that very massive 
GCs contribute to more than 50% to the total 
cluster mass, in fine agreement with the observa- 
tional data of Harris et al. (2006) that indicate 
how up a full 40% of the total mass that is now 
in the GCSs of brightest cluster galaxies is con- 
tributed by massive (present day mass > 1.5 x 10 6 
M Q ). It may be also worth remembered the spe- 
cific paper by Baumgardt et al. (2003) devoted 
to the modeling of Gl, the massive GC in M31, 
obtaining for it a mass 8 x 10 6 M . 

Consequently, one can just make a hypothe- 
sis on the initial abundance of massive clusters 
and see whether it gives results consistent with 
available observations, checking, of course the de- 
pendence of results on the assumption. This was 
done in Capuzzo Dolcetta & Vicari (2005), that 
shows (see their Fig. 11) how few tens of massive 
(M > 10 7 M ) clusters suffice to give an accretion 
rate onto a central galactic black hole high enough 
(few M yr _1 ) to sustain an AGN activity. Given 
all this, our main concern is to check whether or 
not the tidal distortion suffered by the GCs due to 
the halo and the bulge is destructive on the time- 
scale needed by df to dissipate the cluster orbital 
energy. The detailed numerical studies performed 
by Charlton & Laguna (1995); Nordquist et al. 
(1999) give encouraging support to this claiming. 
In this paper, we show the results of detailed N- 
body simulations that take into account df and 
two other concurrent effects: i) the tidal interac- 
tion of a GC with the overall galactic field and ii) 
the collision with another passing- by GC. We con- 
sider quasi-radial orbits for GCs, such to maximise 
the tidal effects. 

As a 'by-product' of this work, tidal tail forma- 
tion and morphology is analysed and a particular 
attention is focused on the presence of overdensi- 
ties (clumps) in the tails. Similar sub-structures 
have been detected in the two tidal tails of the 
globular cluster Palomar 5 (Odcnkirchen et al. 
2001, 2003), while the presence of tails surround- 
ing many other GCs is strongly suggested by var- 
ious observations (Lehmann & Scholz 1997; Testa 
et al. 2000; Leon et al. 2000; Siegel et al. 2001; Lee 
et al. 2003). NGC 6254 and Palomar 12 seem to 
show similar overdensities in their tails, too (Leon 
et al. 2000). Various authors attribute the forma- 
tion of clumps in the tails to strong gravitational 
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shocks suffered by the cluster (Combes et al. 1999) 
and, in particular, to the tidal effect due to com- 
pact galactic sub-structures (Dehnen et al. 2004), 
like molecular clouds or spiral arms. Nevertheless, 
the mechanisms of their origin and formation is 
still unclear and deserve an adequate interpreta- 
tion. More recently, Capuzzo Dolcetta, Di Matteo 
& Miocchi (2005) (hereafter CDM) showed that 
clumpy tails emerge also in GCs moving in a reg- 
ular 'smooth' galactic environment on both quasi- 
circular and more eccentric orbits, suggesting that 
such clumpy sub-structures are related to the lo- 
cal decrease of the stars velocity along the tails 
(see also Di Matteo, Capuzzo Dolcetta & Miocchi, 
2005). 

This paper is organized as follows: in Sect. 2 
the numerical modeling is described, the results 
are shown and discussed in Sect. 3 and conclusions 
drawn in Sect. 4. 

2. The model 

We consider GCs as iV-body systems moving 
within a triaxial galaxy represented by an ana- 
lytical potential (see next Section), including the 
deceleration due to dynamical friction caused by 
the field stars. We studied the evolution of a pair 
of GCs whose center-of-mass (CM) initial condi- 
tions were assigned such to put them along very 
elongated orbits. They start moving on oppo- 
site sides respect to the galactic center and os- 
cillate quasi-radially around it, crossing and col- 
liding each other in the galactic center neighbour- 
hood (see, e.g., the upper panel in Fig. 3). Both 
the choice of very elongated orbits and the pres- 
ence of two clusters were motivated by the wish to 
study the combined effects of the tidal interaction 
and of the collisions between GCs. Indeed, when 
such systems undergo the final stages of orbital 
decay, GC-GC merging and close interactions are 
expected to occur. 

We dealt with four different types of clusters, 
named (a) , (b) , (c) and (d) at increasing order 
of concentration (see Table 1 and Sect. 2. 2). Their 
dynamics was modelled in simulation A, that in- 
volved the pair of clusters (a) and (b) , and in sim- 
ulation B that regarded the more compact clusters 
(c) and (d) . We also considered a further case 
(simulation C) where the dynamics of cluster (a) 
, with the same initial conditions as in simulation 



A, was followed without the presence of the other 
cluster. This latter simulation lasted less than the 
others due to CPU-time limitations. Note that the 
total CPU-time spent by all the simulations pre- 
sented here is about 30,000 hours divided among 
the 64 processors used on an IBM SP4 system 
(granted by the INAF-CINECA agreement). 

In this paper, unless otherwise specified, lengths, 
masses and time are measured, respectively, in 
unit of the galactic core radius n,, of the galactic 
core mass Mb and of the galactic core crossing time 
t b = {rl/GMb) 1 / 2 . Note that the numerical values 
in physical units of the mentioned parameters are 
irrelevant for the results of the simulations, that 
can be always scaled as long as the ratios t/tb, 
r/rb, MjMb (with M the total cluster mass) are 
kept unchanged. This is strictly true only if one 
can neglect, like in our case, the effects due to 
2-body stellar collisions (as discussed in Sect. 2. 2), 
otherwise the dynamics would depend also on the 
mass of the single star in the clusters. 

The reference frame was fixed with the origin 
at the galactic center and the x and z axes, re- 
spectively, along the maximum and minimum axis 
of the triaxial ellipsoid. The cluster CM was ini- 
tially located on the x-axis: cluster (a) and (c) at 
x = —4.15, and cluster (b) and (d) at x = 3.95. 
In all the cases the initial velocity components 
were (0,0.05,0.025). Because of the CPU time 
limitations, we were not able to consider a wide 
set of initial conditions, thus we decided to choose 
initial velocity and position in such a way as to 
set upper limits to the disruptive effects exerted 
on massive clusters by the tidal interaction with 
the galactic field and by the close interaction with 
other GCs. For this reason, we chose almost ra- 
dial orbits so that the GCs cross twice per pe- 
riod the galactic center, where the tidal interaction 
is strongest, and they undergo head-on collisions 
with each other while crossing the center. 

Moreover, the range of initial orbital parame- 
ters compatible with our assumed starting condi- 
tions obviously depends on the actual age of sim- 
ulated clusters, that cannot be univoquely defined 
(clusters of different ages and different initial con- 
ditions may have the same apocenter and veloc- 
ity we chose to start simulations). However, we 
checked that in the adopted triaxial potential all 
GCs with masses in the range (1.5-2.0) x 10 7 M 
(see Section 2.2 and Table 1) and moving on cither 
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box orbits or loop orbits with pericenter ~ and 
initial apocenter in the interval 3.5-4 kpc decay 
(due to dynamical friction) to our starting condi- 
tions in about 1 Gyr (having set Mb = 3 x 10 9 M 
and rt, = 200 pc). The decay time can be easily 
scaled to any GC mass, due to its inverse linear 
proportionality to it, that means a decay time 2 
Gyr for a GC half that massive. 

As regards the computational techniques, we 
adopted a direct A-body representations for the 
stars in the clusters, simulating their dynamics by 
mean of the parallel MPI 'ATD' code (Miocchi & 
Capuzzo Dolcctta 2002) whose main features are 
resumed in Appendix B. 

2.1. The galactic model 

We considered the same self-consistent triaxial 
model described in de Zeeuw & Merritt (1983) and 
in PCV; the same model was also used in CDM 
to study tidal tails formation around clusters in 
absence of df and on not very elongated orbits. It 
corresponds to a non-rotating ellipsoidal, triaxial 
distribution of matter with axial ratios 2:1.25:1 
(Schwarzschild, 1979), which leads to a projected 
profile in agreement with that observed in spirals 
spheroids (see, e.g., Bertola ct al. 1991; Matthews 
& de Grijs 2004) and in elliptical galaxies (see, 
e.g., Wagner 1988; Davies et al. 2001; Statler et al. 
2004). The potential produced can be expressed 
as the sum of a spherically symmetric term due to 
a density following the modified Hubble's law 



Pb(r) = Pbo 
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with r& the core radius and pbo = Mb/rf, plus 
other two non-spherical terms that give the triax- 
ial behaviour, i.e. 
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with A = 4irGM b and 



(3) 

(4) 
(5) 




$ 2 



3z 2 - r 2 



Cl 



2(r 2 + c 2 r 2 ) 3 / 2 ' 



-3c 3 



x 2 - y 2 



(r 2 +c 4 r 2 ) 3 / 2 ' 



The coefficients Ci have been chosen in such a way 
to have density axial ratios roughly constant with 
r (de Zeeuw & Merritt 1983), i.e.: ci = 0.06408, 
c 2 = 0.65456, c 3 = 0.01533, c 4 = 0.48067. We 
refer to the parameter Mb as the galaxy mass, 
though this model yields a total infinite mass. Ac- 
tually, M b ~ 0.45M(rb) being M(r b ) the mass en- 
closed in the sphere with radius r\>. The force pro- 
duced by the galactic potential is evaluated analyt- 
ically and then added to each particle acceleration 
during the simulations. 

With regard to the df, we used the gener- 
alization of the Chandrasckhar formula (Chan- 
drasekhar 1943) to the triaxial case (PCV), with a 
self-consistent evaluation of the velocity dispersion 
tensor, taking also into account that the GC is an 
extended object (see Appendix A for details). We 
checked the dependence of df decay on applying it 
to the center of mass of the GC and to its center 
of density (see Sect. 3.2). 

2.2. The cluster model 

Our A-body simulations involved four different 
clusters. They sample a King distribution, with 
M the total mass, a the central velocity disper- 
sion, rt and r c , the limiting and the King radius 
respectively, c = log(r t /r c ) the concentration pa- 
rameter and, finally, t c = (r^/GM) 1 / 2 the core- 
crossing time. The 'limiting radius' is the radius 
at which the King distribution function drops to 
zero to model the presence of the external field 
(King 1966). Of course, for consistency, the lim- 
iting radius should be of the same order of the 
tidal radius corresponding to the local tidal field. 
Anyway, in order to model the presence of some 
tidal debris around the clusters since the initial 
time, the limiting radius was chosen to be 20-60% 
larger than the maximum tidal radius. This max- 
imum value is achieved at the apocenter, i.e. at 
the initial cluster position, where it can be esti- 
mated as r (M/M(r )) 1/3 ~ 0.3 (being r ~ 4 
and M(r ) ~ 14). However, the amount of cluster 
mass lying outside the tidal radius is in any case 
less than about 1% of the total mass. The initial 
values of the parameters are listed in Table 1. 

Each cluster was represented with N = 5 x 10 5 
'particles' in the numerical model. The particles 
were assigned a mass according to a Salpeter's 
mass distribution (dN/dm cx m -2 ' 35 ) in the range 
3.3 x 10" 11 + 3.3 x 10~ 9 (if, e.g., M b = 3 x 10 9 M 



5 



then this range is 0.1 -r 10 M Q ) then, all the par- 
ticle masses were uniformly re-scaled such to give 
the desired M. The choice to represent the GC 
with the 'right' value for the total mass, thus over- 
estimating the individual star masses for a factor 
Ntme/N, where N trU e is the real number of stars 
in the GC, is the same done by, e.g., Combes et 
al. (1999). Anyway, both the relative stellar abun- 
dances and the total binding energy are correctly 
reproduced, thus we expect that the representa- 
tion of the external tidal effects on the internal ve- 
locity distribution (which actually drives mass loss 
from the cluster) is correct as long as the spurious 
heating effect introduced by the overestimate of 
stellar masses is not relevant. As a matter of fact, 
we checked (see later in this Sect, and Sect. 3.1.3) 
that this heating is absolutely negligible over the 
time length of our simulations. 
There are, of course, other possibilitcs to pick par- 
ticle masses. For example, Ideta & Makino (2004) 
and Dehnen et al. (2004) assume all particles to 
have the same mass. This corresponds to an in- 
dividual mass of about 130 M Q in the case of the 
simulations of the formation of the massive ui Ccn 
cluster (Ideta & Makino 2004), and of about 1 
M Q in the case of the study of tidal tail formation 
around the light Pal 5 cluster. Another possible 
choice is that of having particle masses distributed 
according to an assumed stellar mass function, 
subsequently rescaling other relevant GC charac- 
teristic parameters (as done by, e.g., Baumgardt 
& Makino (2003)). 

All these choices have critical aspects, unavoidably 
clue to the limited number of particles allowed. In 
particular, as said above, our choice may imply 
spurious heating effects (as shown in the context of 
numerical galaxy formation by Stcinmctz & White 
(1997)), in what fluctuations over the mean field 
are enhanced by the too high 'star' mass value. 
However, fluctuations should not be exceedingly 
important because the total mass of the system 
(the most relevant parameter in determining the 
GC mean field) is kept at its right value. In this 
respect, note that the half-mass relaxation time 
t r h can be evaluated as (Spitzer 1987): 



where, e.g. for the most compact cluster (d) , 
t c /t b ~ 3.7 x 10~ 2 (Tabic 1). Following Giersz & 



Heggie (1994), A ~ 0.1A, thus the above written 
formula gives t r h ~ 240ib for the system (d) in the 
numerical model. Since our simulations reached 
about AOtb, we can say that 'spurious' collisional 
effects should be always negligible. We say 'spuri- 
ous' because, of course, such effects are even more 
negligible for the 'real', physical, clusters. Indeed, 
if, say, Mb = 3 x 10 9 M Q , then cluster (a) , for 
instance, would have M = 2 x 10 7 M Q ; hence, 
if the stellar average mass is <m>~ 0.3 M©, it 
would be made up of N = 7 x 10 7 stars, making 
the real relaxation time even longer than the sim- 
ulated time. Finally, our simulations start with 
clusters that have an age less than their internal 
2-body relaxation time, so we assumed no initial 
mass segregation. 

In Fig. 1 the orbits of the most concentrated 
clusters, (c) and (d) arc plotted, as turned out 
from simulation B (for which clusters keep always 
a well defined core). They are quite elongated 
box orbits and represent the motion of the clus- 
ter center- of- density (CD), i.e. the average of the 
particles positions weighted with the local density 
instead of the mass (Casertano & Hut 1985). In 
most cases, we decided to take the CD as the best 
suitable reference for that regards the study of the 
internal cluster properties. Indeed, as we verified, 
the CM position is strongly influenced by the very 
extended tidal tails which quickly formed, so to be 
located even well outside the clusters' core. 

3. Results 

3.1. The effects of the tidal shock on the 
GCs internal structure 

3.1.1. The clusters internal 'heating' 

In Fig.s 2 and 3, for each cluster we plot the 
total internal kinetic energy (K) and the internal 
gravitational potential energy (U) of that part of 
the cluster that is enclosed in the sphere centered 
at the CD with a radius equal to the initial King 
radius of the cluster itself. 'Internal' means that 
K is evaluated with respect to the cluster CM, 
while in U the 'external' potential produced by the 
galaxy and by the other cluster is not taken into 
account. Considering the least compact clusters, 
(a) and (b) , one can note, from Fig. 2, that the 
potential energy drops off (whereas K shows an 
'impulse') when the clusters pass across the core 
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of the galaxy, in agreement with the fact that their 
inner part undergoes a violent compression at that 
moment, as clearly confirmed by Fig. 4. After- 
wards, the potential energy seems to recover its 
previous level, but not completely because after 
each core-crossing a certain amount of potential 
energy is irreversibly lost by the clusters, in favour 
of some internal 'heating'. All this is true espe- 
cially for the least compact clusters (a) and (b) , 
while (c) and (d) models suffer this tidal influence 
in their external regions only. Note, indeed, the 
much smaller variations in U and K (Fig. 3) and 
the nearly constant Lagrangian radius of 10% of 
the mass (Fig. 4) for the most concentrated clus- 
ters. 

The energy input is due to the tidal interac- 
tion with the environment and not to the direct 
collision between the two clusters, as verified by a 
careful comparison of the time behaviour of the x 
coordinate of the CDs and of the energy (Fig.s 2 
and 3). This is confirmed by the results of the sim- 
ulation C that involves the presence of the cluster 
(a) alone, thus excluding any effect due to GC-GC 
collision. It can be affirmed that, at least within 
the simulated time, both the energy behaviour and 
the motion of the cluster show no appreciable dif- 
ference when it is alone in comparison with the 
case when it is in a pair. 

Such a conclusion is not that surprising for or- 
bits, like those investigated here, that give rise to 
collisions with a high relative velocity V (V/a <~ 
60, being a ~ 0.1 the cluster internal velocity 
dispersion). Notice, indeed, that the duration of 
a head-on collision is ~ r c /V <~ 10 ~ 2 which is 
shorter than the crossing time of our simulated 
clusters, so that the 'impulse approximation' can 
be applied in order to estimate the change in the 
clusters energy (per unit mass), E = U+K, due to 
one collision, that is (Binney & Tremaine 1987): 
AE ~ G 2 M 2 {V 2 r 2 c )- 1 ~ a A V~ 2 ~ \E\(a/V) 2 . 
Hence the number of collisions needed to change 
appreciably the internal energy of the clusters is 
of order of \E\/AE ~ (V/a) 2 ~ 10 3 , i.e. much 
more than that occurred in our simulated time. 
Of course, for smaller V such an approximation 
is no longer valid and the close encounters may 
contribute furtherly to enhance the probability of 
merging between clusters, as suggested by the sim- 
ulations of superclustcrs merging in Fellhauer & 
Kroupa (2002) and discussed in detail in a forth- 



coming paper. 

Note that while the loosest clusters (a) and (b) 
appear almost destroyed at t > 10, when E > 
and the CD oscillation around the galactic center 
begins to lose its regularity (see Fig. 2), the more 
compact clusters in the simulation B are able to 
survive the tidal stress (Fig. 3). 

3.1.2. The mass loss 

The mass loss of the simulated GCs is due to 
the tidal interaction with the galactic potential. 
Indeed, we checked that the same cluster mod- 
els used in the simulations if evolved at rest and 
without any external field for a time equal to the 
time of the simulations, exhibit a mass loss rate 
(in terms of particles escaping out of the limiting 
radius) not greater than 1%, presumably due to 
the (small) relaxation effects. 

We quantify the amount of the mass lost adopt- 
ing various criteria. According to an 'energetic' 
one, /iE is the fraction of the mass lost made up 
of the stars whose internal energy — i.e. that due 
only to the other stars of the same GC and evalu- 
ated in the reference frame with the origin at the 
cluster CD — is non-negative. This criterion is not 
rigorous because the sign of the individual star en- 
ergy does not guarantee its 'dynamical destiny'. 
Another indication of the mass lost from the clus- 
ters may be given with reference to observations. 
Defining as 'lost' those stars that belong to re- 
gions where the local GC density falls below apt,, 
the choices of a = 1 and a = 0.1 lead to the quan- 
tities hi and Ho.i as fractions (to the total initial 
"observable" GC mass) of mass lost from the GC 
reported in Fig. 5. 

Note that all the escaping criteria give very sim- 
ilar behaviours, especially for the most compact 
clusters. This behaviour is well fitted by a law 
of the form fj, <~ 1 — acxp(— t/r); in Table 2 the 
best-fit parameters are given for the various GCs 
according to the different criteria adopted. Note, 
also, that although the initial limiting radii are 
greater than the local tidal radius, the cluster mass 
initially outside the tidal radius is, in any model, 
less than the 1% of the total mass, giving a negli- 
gible contribution to the overall mass loss. 

Even if the statistics is poor (only 4 mod- 
els), nevertheless we quantified a correlation be- 
tween the King concentration parameter c and 



7 



te (the Pearson's linear correlation coefficient be- 
tween lnT£ and lnc is 0.992). The best power law 
correlation in the bi-logarithmic plane (least % 2 fit, 
corresponding to x 2 = 0.091) is teI% ~ 14 x c 61 . 

3.1.3. Density profiles and mass segregation 

For the two most compact clusters ((c) and (d) 
) in the simulation B, we were able to fit the in- 
ner radial density profile with a King distribution 
throughout the whole simulated time, as shown in 
Fig.s 6-7. We also studied the time behaviour of 
the concentration parameter (c) , of the central ve- 
locity dispersion (<r) and of the limiting (r t ) and 
King (r c ) radii for both clusters (Fig. 8). Four 
points are of particular interest: i) the GCs remain 
bound and keep a well defined spherical "core" ; ii) 
the inner part of the clusters tends towards less 
and less concentrations (see also the bottom panel 
of Fig. 5); iii) c shows wide oscillations coincid- 
ing, temporally, with the repeated compressions 
and re-expansions experienced by the clusters at 
each passage in the galactic core (see Fig. 4); iv) 
a appears to be roughly constant. The decrease 
of the concentration is due to both the growth of 
the King radius and to the decrease of the limiting 
radius. The former is caused by the continuous re- 
virialization of the core as the cluster loses mass, 
while the latter is a direct consequence of the tidal 
erosion occurring preferentially in the system out- 
skirts. 

It is worth noting that in our case the tidal ero- 
sion is strong enough to act much more rapidly 
than 2-body collisions (r < t r h/10). This ex- 
plains the apparent contradiction of our results, 
in particular point (ii), with those of other au- 
thors investigating the role of tidal interaction, 
e.g. Gnedin et al. (1999). Indeed, they found, 
by means of Fokkcr-Planck simulations of a GC 
undergoing tidal forces much weaker than in our 
case, that the tidal erosion acts just to acceler- 
ate the core collapse on its characteristic time 
scale, which is of the order of tens of t r h- At this 
same regard, we note that an initial growth of the 
core of not too concentrated clusters was already 
found by Spitzer & Chevalier (1973) and Spitzer & 
Shull (1975) (in the cases of single-mass and mul- 
timass Montecarlo models, respectively) on a time 
scale short compared to t r h, followed by a rapid 
evaporation-induced core collapse. 

In order to investigate possible mass segrega- 



tion phenomena — not expected to be relevant be- 
cause the simulated time is short in comparison 
with the 2-body relaxation time — we analysed the 
evolution of the average mass of stars within three 
different clusters regions in simulation B: only for 
cluster (c) in the pair, a rather small increase of 
the average mass is found (~ 10% of the initial 
value) in the inner part, i.e. in the sphere where 
the cluster included initially 20% of its mass. In 
the outskirts and in the tails no appreciable mass 
segregation occurs. This confirms the absence of 
significant collisional effects that, if present, would 
have indicated either an influence of the external 
tidal field in decreasing the relaxation time, or a 
too small number of particles with respect to the 
actual number of stars. 

3.2. The orbital decay 

In the framework of our main scientific moti- 
vation, it is important to quantify how fast GCs 
decay towards the center of their parent galaxy. At 
this purpose we defined the adimensional quantity 

, v E orb (t) - W 
€orb(t) = — s -r , (7) 

£,0 — W 

where E orb is the orbital energy of the CD of the 
cluster, E = E orh (0) and W = -4:irGM b /r b is 
the galactic potential well. If neither dissipation 
nor 'energy injection' occur, then £ or f, = 1. In 
Fig. 9 we give the time behaviour of £ or b for the 
clusters (c) and (d) , only, because they keep a 
rather defined core on the whole duration of the 
simulation. In order to understand whether the 
tidal interaction with the galactic field influences 
appreciably the orbital decaying, it is worth com- 
paring the clusters £ or b{t) evolution with that of 
two point-masses, moving within the same galac- 
tic potential with the same initial conditions of the 
clusters and experiencing the same df (Fig. 9), but 
free from tidal effects. Note that in the evaluation 
of the orbital energy of a cluster, the potential en- 
ergy given by the interaction with the other is not 
taken into account; as we said in Sect. 3.1.1, the 
overall evolution of a cluster is not influenced ap- 
preciably by the presence of the other GC. Thus, 
to simplify the comparison, we considered the two 
point-masses as non-mutually interacting objects. 
Otherwise, we would have had to smooth the mu- 
tual gravitational force with a suitably variable 
smoothing radius, in order to simulate the variable 
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tidal distortion occurring for the clusters. More- 
over, the mass of the point-masses is kept constant 
for consistency with the choice of constant clusters 
half-mass as a parameter in the df formula (see 
Appendix A). 

We expected a different behaviour between the 
two cases, mainly because of the tidal interac- 
tion with the galaxy, that occurs only in the full 
A-body case and that, presumably, gives rise to 
a partial "thermalisation" of the orbital energy 
among the internal degrees of freedom, resulting 
in a further form of dissipation. This prediction 
is confirmed by the energy evolution, in the sense 
that the rate of dissipation is greater in the case 
of extended objects than for the two point-masses. 
The peaks shown by £ or f, are due to the strong ac- 
celeration of the clusters CM during their close en- 
counters (the mutual potential energy is not taken 
into account). To quantify the further frictional 
effect on the extended bodies, we followed the de- 
cay of the point-masses up to the time tdec when 
£ or & ~ 10~ 2 , i.e. when their orbit is confined to 
a region of size comparable with that of a typi- 
cal GC. This happens at tdec ~ 400 for both the 
point-masses. On the other side, from Fig. 9, we 
can extrapolate a value of tdec f° r the A^-body sys- 
tems assuming a constant energy decay rate, as 
tdec 5; 180 for the cluster (c) and td ec ^ 280 for 
the more compact but less massive cluster (d) . 
Actually, these are overestimates, because the en- 
ergy decay rate is generally increasing with time 
(see PCV and Capuzzo Dolcetta & Vicari 2005). 

We must note, however, that the energy dis- 
sipation of the clusters depends on the particu- 
lar way the df is evaluated. Indeed, as described 
in Appendix A, we used, in the Chandrasekhar 
formula, the kincmatical quantities of the clusters 
CM that, due to the quick formation of large tails 
of stripped material, exibits a very rapid decay to- 
wards the center, where the density and velocity 
dispersion of the galactic model have larger val- 
ues than at the location of the main body (core) 
of the cluster. Actually, without an A-body self- 
consistent representation of the galaxy in which 
the satellite moves, the way to compute and as- 
sign the df deceleration to the bodies of a very 
extended object is not a trivial issue. It is diffi- 
cult to model correctly the df changes induced by 
the cluster distorsion and, maybe, even more dif- 
ficult to take into account the gravitational feed- 



back of the cluster on the galactic nuclear region, 
thus we were forced to do unavoidable simplifica- 
tions. As explained in Appendix A, we assumed 
that the 'global' effect of the df is uniformly dis- 
tributed to every GC star, and we evaluated it as 
if the cluster were concentrated in the CM posi- 
tion with the CM velocity. Even if questionable, 
this is a logical choice that would deserve a dis- 
cussion, which is out of the purposes of this pa- 
per. In any case, we re-simulated (though with 
A = 10 4 , for obvious computational convenience) 
the GCs evolution for models (c) and (d) (as in- 
dividual systems in the galactic field), and with 
df computed by mean of the formula given by Eq. 
(Al) "centered" in the CD instead of the CM, and 
the comparison is shown in Fig. 9. The £ or f, be- 
haviour in these new simulations is flatter than 
that of the case with df evaluated at CM. The 
smaller effect of df is explained by that the CM 
is systematically closer to the galaxy center than 
the CD. Nevertheless the different amount of df 
does not modify qualitatively the global GC evo- 
lution at least within the simulated time. Note 
that a reliable estimate of the mass loss in the 
df-on-CD case cannot be achieved because the re- 
duced number of particles (A = 10 4 ) makes the 
2-body relaxation time even shorter than the sim- 
ulated time (see Eq. 6), thus affecting significantly 
the evolution with (spurious) collisional effects. 

Considering the limited computational re- 
sources available, we chose to adopt an as more 
accurate as possible A-body representation of the 
cluster to study in detail the tidal disruption pro- 
cess, even if this forced us to employ an analytic 
single-component model for the galaxy, and so an 
analytic treatment of the df effect. However, the 
presence of velocity anisotropy for the field stars 
has to be take into account, since it is important 
in altering the efficiency of the df as was proved by 
Binncy (1977) in the axisymmetric case, then by 
PCV in the triaxial case and recently confirmed 
by the numerical simulations of Peharrubia et al. 
(2004). Thus, we decided to use the PCV gener- 
alisation of the Chandrasekhar formula, with the 
self-consistent implementation of the analytic ex- 
pression for the velocity dispersion tensor provided 
by the galactic model we adopted. Moreover, we 
took into account the change in the df caused by 
the non- uniform density of the field stars, through 
a suitable local estimate of the Coulomb log (see 
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Appendix A). 

3.3. Tidal tails morphology 

Tidal tails rapidly form around the clusters and 
follow closely their orbit (see Fig. 10). This tail- 
orbit aligncment has been recently observed for 
the GC Palomar 5 (Odenkirchen ct al. 2001, 2003) 
and also reproduced in various simulations (Charl- 
ton & Laguna 1995; Combes et al. 1999; Dehnen et 
al. 2004; Capuzzo Dolcetta et al. 2005). Since the 
velocity dispersion in the cluster outskirts (where 
most of the evaporated stars come from) is much 
lower than the cluster orbital velocity it is not 
surprising that, there, stars move on orbits simi- 
lar to that of the cluster itself, as clearly shown by 
Fig. 10. It can be also seen a significant spread of 
escaped stars around the apocenter of the cluster 
orbit. This because of the small differences among 
the velocities of the stars at the moment they leave 
the cluster become unimportant at the apocenter, 
where the orbital velocity is very low. 

3.3.1. Clumps and ripples formation 

In our case, where GCs orbits are quasi-radial, 
centrifugal and Coriolis' forces are not important 
in determining the tails shape distorsion investi- 
gated in CDM on a sample of less elongated or- 
bits. Nevertheless, we do observe in our simula- 
tions the presence of stellar overdensities along the 
tidal tails, similar to those seen around the men- 
tioned Pal 5 and in the simulations in CDM. 

The cluster (b) is plotted at various times close 
to the first passage at the rightmost apocenter, 
in Fig.s 11-12. A "ripple" starts to form around 
x = 3.1 at t = 8.6, as a 'wave-like' overdensity. 
Then, the cluster travels for a while across and 
then above the formed ripple (t = 9.2-9.7), before 
reaching its apocenter. It can be also seen that 
at t = 10.3 (Fig. 12) another ripple forms again 
around x = 4.1, while at t — 9.7 the one previ- 
ously formed begins to move inward giving rise, 
at t = 10.5, to a clump — an overdensity with a 
roughly spherical shape — around x — 3.2. All 
these plots are projections on the xy plane, how- 
ever all the structures are almost aligned along 
the orbit around which (i.e. around the x-axis) 
they show a roughly cylindrical symmetry; as an 
example the ripples are, actually, thin discoidal 
structures. This can be seen in the available 3-d 



animations (see below). 

To explain the overdensities found in the tails 
of the mentioned Palomar 5, Dehnen et al. (2004) 
raised an hypothesis stating that clumps could be 
due to the effect of the interaction with Galac- 
tic substructures (like giant molecular clouds, spi- 
ral arms, dark-matter sub-halos or massive com- 
pact halo objects). Basically their opinion is that 
small-scale overdensities like clumps can only be 
built up by the tidal interaction with fluctuations 
on relatively small-scale of the external gravita- 
tional field. On the contrary, the results shown 
here seem to confirm the findings of CDM, in the 
sense that in the absence of any small-scale sub- 
structure in the overall potential clumps form too, 
thus suggesting that a further mechanism is at 
work. 

Analysing the local velocity measured along the 
orbital path and the energy of the stars in the clus- 
ter field, we can affirm that both clumps and the 
apocenter ripple are not gravitationally-bound ag- 
gregates but, rather, overdensities due to the local 
deceleration of the stellar motion, as already ver- 
ified in CDM and in Di Matteo et al. (2005). On 
the contrary, where the velocity increases towards 
the direction of the stellar mean motion along the 
tails, the tails tend to rarefy and an underdensity 
occurs. The metaphor of the "motorway traffic 
jam" mechanism used in Dehnen et al. (2004) is 
effective here, as can be seen in Fig. 13. The cause 
of the overdensity is the strong deceleration of the 
stellar flux immediately on the left of the ripple 
location (around x = 3.25). Note that inside the 
cluster the profile of <v x > is nearly flat as ex- 
pected for a bound (almost rigid) object, while in 
the ripple region it exhibits a non-zero gradient. 

It is interesting to compare a configuration like 
those shown in Fig. 12 at t — 10.3 and t — 10.5, 
with the 'arcs' of material observed around the 
spectacular shell galaxy NGC 3923 (Fort et al. 
1986; Pence 1986) as well as with the results 
of the simulations in Hcrnquist & Quinn (1987, 
1988); Heisler & White (1990) and Bournaud ct al. 
(2004). Although different scenarios and length- 
scale are involved, the qualitative features of the 
debris produced by the tidal interaction of a com- 
pact and smaller object with a larger density dis- 
tribution, are rather similar to what we found in 
this paper. In particular, sharp-edged structures 
in form of 'ripples' or 'arcs' are formed when the 
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satellite moves on a very elongated orbit in a reg- 
ular potential generated by a much more massive 
and larger object. 

To conclude this Section, we give the web 
address http : //inaf . cineca. it/KP/kp4 . html 
where the reader can get animations regarding 
some parts of the simulation A. The first anima- 
tion refers to the whole simulation and one can see 
how the least compact GCs get destroyed. The 
second animation is a zoom around the apogalac- 
ticon, that corresponds to the Figures 11-12 of 
the sub-structures formation. 

4. Conclusions 

In this first of a series of papers we investigated 
the effects of the central galactic environment on 
the structure of globular clusters decayed in the 
inner region because of dynamical friction. Our 
main concern is to answer the following questions: 
i) to what extent are clusters able to survive the 
strong tidal interaction with the galaxy? ii) do 
they keep their compact structure long enough to 
permit dynamical friction to dissipate completely 
their orbital energy? iii) what are the effects of 
close encounters and collisions between clusters? 

The main conclusions of our work can be sum- 
marized as follows. 

Sufficiently compact clusters (initial King con- 
centration parameter c > 1.2) can survive the 
tidal interaction with the external fields for, at 
least, the duration of our simulations, i.e. for 8 
passages across the galactic center (i ~ 404, be- 
ing 4 the galaxy core crossing time), while those 
with c = 0.8-0.9 are almost disgregated after only 
2-3 passages (~ 10-154). The mass loss from 
the clusters as a function of time follows an ex- 
ponential law ~ e~ t/,T with t up to ~ 704 for 
the most compact cluster (c = 1.3). Even if still 
without a high statistical reliability, we found a 
correlation between the mass loss time-scale and 
the concentration parameter that can be quanti- 
fied as t ~ 14c 61 . This means that a GC with an 
initial c > 1.6 keeps bound a substantial amount 
of its mass up to t ~ 3004, i-c. up to the complete 
orbital decaying. 

During their evolution, the survived clusters 
maintain the initial King profile in the inner re- 
gion, with a decreasing concentration and a con- 
stant central velocity dispersion. A small degree of 



mass segregation occurs at the end of the simula- 
tions: the average stellar mass in the central clus- 
ters region increases of about 10% of the initial 
value, but the fluctuations are even larger. The 
very low level of mass segregation is not surpris- 
ing, for the half-mass relaxation time is ~ 6 times 
longer than the simulated time. 

With regard to the orbital decay, we found that 
the tidal interaction with the field gives rise to a 
dissipation of the orbital energy with a rate com- 
parable to that given by the df, when this is eval- 
uated at the cluster CM. However, even when the 
CD is used instead of this latter, we found that 
the tidal interaction provides a further important 
mechanism of orbital decaying besides dynamical 
friction (included in our simulations). This is im- 
portant for the validity of the nucleus accretion 
model according to which the nuclear region re- 
ceives a large amount of mass in form of clusters 
that have lost their orbital energy (Capuzzo Dol- 
cetta 1993; Capuzzo Dolcetta & Vicari 2005). 

Another relevant finding of this work is that, 
in the cases considered here, the tidal interaction 
between the two clusters, even during face-on col- 
lisions, produces negligible effects on both their 
internal evolution and their orbital evolution (at 
least for t < 404), confirming the estimates done 
according to the impulse approximation. 

All these results indicate that sufficiently mas- 
sive and compact clusters can survive the disrup- 
tive effect due to the strong tidal forces exerted 
by the central environment. Dynamical friction 
and tidal dissipation are then able to brake many 
globular clusters before they are disgregated, so to 
allow the formation of a dense and massive super- 
cluster resulting from merging events among clus- 
ters. The simulations shown here did not achieve 
the final merging stage between the objects, be- 
cause of the limits in the computational resources. 
In a forthcoming paper, we will present results of 
a new series of simulations directly regarding the 
merging among clusters. Indeed, the importance 
of the present work is also to get physically reli- 
able and realistic initial conditions from which the 
simulations of the final merging process will start. 

A further result of this work concerns the for- 
mation and structure of tidal tails. Tidal tails 
form rather quickly, with a clear tendency to align 
along the orbital path. They exhibit complex mor- 
phology, with the presence of clumps and 'rip- 
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pies' as those observed around globular clusters 
and shell galaxies, too. Such overdensities are 
not gravitationally bound aggregates but, rather, 
their formation seems to have a kinematical origin, 
being connected to the deceleration of the stellar 
'flux' motion along the tails. 
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A. Treatment of the dynamical friction 

The effect of braking due to df on globular clusters was followed by the same approach used in PCV. The 
classic formula developed by Chandrasekhar (1943) for the df deceleration term has been extended by PCV 
to the triaxial case, in partial analogy with the (Binney 1977) extension to the axisymmetric case, obtaining: 

&df = -7iViei - 72^2 - 73^363 (Al) 

where ej (i = 1, 2, 3) are the eigenvectors of the velocity dispersion tensor of the galaxy stars and Vi is the 
component of the velocity of the baricenter of the GC along the e, axis. The coefficients 7$ are (see PCV): 

7l = ^W^MlnA x (A2) 




where: p(r) is the mass density of background stars, In A is the Coulomb's logarithm, M is the mass of 
the test object, G is the gravitational constant, Ui is the eigenvalue, corresponding to e,, of the velocity 
dispersion tensor evaluated in r, and u is the ratio between Ui and the greatest eigenvalue, set as o\. 

The velocity dispersion was computed and presented by Merritt (1980) for the Schwarzschild ellipsoid, 
both in the case of a rotating and a non-rotating model. In this paper we considered just the non-rotating 
model. 

For computational convenience, the Merritt's data were fitted obtaining three analytical expressions for 
the eigenvalues cr 3 < (T 2 < cri at f = |r|/r h : 

a\ = 3.1e- f / 9 - 2 

° l = l + 0.43fi-60 (A3) 
2 



1 + 0.44f 



=1.7 



(expressed in unit of GM^jr^). Analogously, the following fitting formulas for the Euler angles giving the 
orientation of the local reference frame where is diagonal were determined: 

2f 2 ' 3 

a = 7T / 2+ 1 +02 f2.3 [arccos(z/f) - tt/2] 

n 2f 2 - 25 ( v \ 

^"/ 2 + TTa*-" csi °(^Fj <A4 » 

0.8f 3 ( ij \ 

7 = ; , no -a arcsin 



l + 0.8f 3 V V* 2 +y 2 J 

The fitting formulas A3 and A4 are a slight improvement of those reported in PCV. 

Given 6m and b m respectively, the maximum and the minimum impact parameter of the test object 
with field stars, we adopted a variable and local Coulomb's logarithm In A = \n(bM /b m ) assuming b rn = 
GM 1 / 2 /(J 2 (r), where cr 2 (r) = erf + a\ + a\ and M 1 / 2 — M/2 is the cluster initial half-mass. Indeed, since 
the GC is not a point-mass, the df is made less efficient by the 'weakening' of the encounters with impact 
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parameters smaller than the size of the object itself, assumed as the half-mass radius (see e.g. Binney & 
Tremaine 1987). Of course, this reasoning is based on the hypothesis that at least the core of the cluster 
survives to the tidal stress so to make df to be able to continue acting on a compact object with mass 
~ M1/2, throughout the whole simulated time. Our choice of letting 6m = lOrj seems quite adequate to 
our quasi-radial orbits, noting that the stellar density at a distance r > lOr^ from the galactic center falls 
below 10~ 3 times the central value. Choices of bu as function of the galactocentric distance, like for instance 
bu = t (Hashimoto ct al. 2003), are of not straightforward use in our cases. 

Finally, the acceleration given by Eq. (Al) is evaluated either at the CM or at the CD position of the 
cluster and then added to the accelerations of every particle of the cluster. 

B. The TV-body approach 

The simulations have been performed by means of the code 'ATD' (Miocchi & Capuzzo Dolcetta 2002). 
It is a a parallel tree-code that follows, in part, the Barnes & Hut (1986) algorithm. For the potential 
due to distant 'particles' it uses a multipolar expansion truncated at the quadrupole moment and has been 
parallelized to run on high performance computers via MPI routines, employing an original parallelization 
approach. A good resolution in the evaluation of the total gravitational force acting on each star of the 
cluster is achieved by a direct summation of the contribution given by neighbours. To avoid instability 
in the time-integration, the 1/r gravitational potential has been smoothed by using a continuous /3-splinc 
function for r < e, such to give an exactly Newtonian potential for r > e (Hcrnquist & Katz 1989) and a force 
that vanishes as r for r — > 0. The time-integration of the 'particles' trajectories is performed according to the 
leap-frog algorithm. This latter uses individual and variable time-steps according to the block-time scheme 
(Aarseth 1985; Hernquist & Katz 1989), in addition with corrections we implemented in order to keep the 
same order of approximation also during the time-step change. Denoting by t c the minimum core-crossing 
time of our simulated clusters, the maximum allowed time-step is Ai max = 0.01£ c , while the minimum is 
Ai m i n = Ai max /2 10 , thus fastest particles may have a time-step as small as <~ 10 _5 t c . To choose the time-step 
of the i-th particle we found that the best compromise between accuracy and speed in the time-integration 
is via the formula: At* = 0.05 x min{ (di/aj) 1 / 2 , di/vi}, where Vi is the velocity of the particle relative to its 
first neighbour (or to the CM of its closest box), di its distance from the first neighbour and Oj the modulus 
of the total acceleration of the i-th. particle itself. 

In all the runs the softening length e is set to 10~ 5 r{, (i.e. 2 x 10~ 3 pc if = 200 pc), so to have 
(e 3 /GM) 1 / 2 <~ Ai m i n . This is a minimal condition to ensure that the level of accuracy in the forces evaluation 
agrees with the accuracy in the time-integration scheme. Note that such an e is much smaller than the typical 
intcrparticlc distance so to keep the overall correct Newtonian behaviour. For N lower than the real number 
of stars, as in our case, the choice of the gravitational smoothing radius is a matter of compromise between 
the need for mantaining the Newtonian behaviour (the smaller is c the more 1 correct is the; force: evaluation) 
and the need for avoiding a spurious increase of collisionality (the larger is e the closer is the model relaxation 
time to the real one). As a matter of fact, the results discussed in Sect. 3.1.3 suggest that e is large enough 
not to introduce fictitious collisional effects. 
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Fig. 1. — Center-of-density trajectories on the xy 
plane, for cluster (c) (dotted line) and (d) (solid 
line) in the Simulation B. 



This 2-column preprint was prepared with the AAS IATfrjX 
macros v5.0. 



17 



Table 1 
Clusters initial parameters. 



cluster model 


M x 10 3 


n 


c 


r c x 10 2 


t c x 10 2 




Simulation 


a 


6.7 


0.47 


0.8 


7 


23 


0.13 


A, C 


b 


5 


0.36 


0.9 


4.5 


13 


0.13 


A 


c 


6.7 


0.49 


1.2 


2.7 


5.4 


0.14 


B 


d 


5 


0.38 


1.3 


1.9 


3.7 


0.14 


B 



Note. — Parameters list for the initial cluster models. Masses, lengths and time are 
in unit of Mb, H and tb respectively. 



Table 2 

Parameters for the mass loss fits. 



cluster model 


as 


TE 


ao.i 


To.l 


ai 


Tl 


a (sim. A) 


3.8 


3.2 


0.75 


35 


0.51 


15 


b 


1.9 


9.7 


0.91 


33 


1.2 


14 


c 


0.77 


46 


0.83 


56 


0.72 


45 


d 


0.89 


65 


0.88 


74 


0.79 


69 



Note. — Parameters of the mass loss evolution 1 — 
ae~ l l T , according to the criteria based on: the total 
internal energy (oe,te), the density contrast of 10% 
(fflo.i)TD.i) an d 100% (ai,n). Time is in unit of tb- 
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Fig. 3. — From top to bottom: CD motion along 
z-axis for the pair of the more compact clusters 
((c) and (d) ) in the simulation B, time behaviour 
of the potential, kinetic and total internal energy 
(solid line: model (d) , dashed: model (c) ). En- 
ergies are evaluated as in Fig. 2. 



Fig. 2. — From top to bottom: CD motion along 
z-axis for the pair of (looser) clusters (a) and (b) 
in the simulation A; time behaviour of the poten- 
tial (U), kinetic (K) and total (E) internal energy 
(solid line: model (b) , dashed: model (a) ). En- 
ergies are evaluated in the center-of-density frame 
and they are concerned only with the stars in the 
sphere with radius equal to the initial King radius 
of each cluster. For the loosest cluster (a) the en- 
ergy is plotted up to t ~ 24 because, afterwards, 
it is almost destroyed. The case of the cluster (a) 
evolving alone in the simulation C showed no ap- 
preciable differences in respect to the simulation 
A. 




20 30 40 



Fig. 4. — Lagrangian radii evolution for the 4 clus- 
ter models (as labelled). In each panel, the radii 
of 10, 50 and 90% of the total mass are shown. 
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Fig. 5. — Different evaluations of the time evolu- 
tion of the fraction of mass lost from the clusters 
(see text). Bottom panel: time evolution of the 
central density normalized to the galaxy central 
density. Solid line: model (d) , dotted: (c) , short 
dashed: (b) , long dashed: cluster (a) in the sim- 
ulation A. The curves are time-averaged over 3tf,. 
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Fig. 6. — Cluster density profile (squares), for the 
model (d) at various times. The sampled spher- 
ical shells have width such to contain a constant 
number of particles. The density is normalized to 
the galaxy central density. The solid curve is the 
best King model fit. 
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Fig. 7. — As in Fig. 6 for the cluster (c) . 
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t 

Fig. 8. — From top to bottom: evolution of the 
concentration parameter, c, of the central velocity 
dispersion, a, and of the King and limiting radii 
(r c , r t ) for the best King model fits to clusters (c) 
(open squares, dotted line) and (d) (solid squares 
and solid line). 



Fig. 9. — Time behaviour of the fraction of or- 
bital energy kept, at time t, by the clusters in the 
simulation B (solid line) and by the corresponding 
point-mass (dotted line). Upper panel: cluster (c) 
; lower panel: cluster (d) . For comparison, the 
case with the df evaluated at the cluster CD by 
mean of the 'reduced' simulations (see text) is also 
plotted (dashed line). 
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Fig. 10. — Configuration of cluster (b) at t = 17. 
The projected trajectory of its center-of-dcnsity 
motion is also plotted. The y-axis scale is ex- 
panded to show better the alignment between the 
cluster orbit and the tails. 




-0.2 r t=9.2 



12 3 4 



Fig. 1 1 . — Cluster (b) is sketched at various times 
(as labelled) before the first passage to the right- 
most apocenter (galactic center is at x = 0). 



Fig. 12. — Continuation of Fig. 11 for the same 
cluster (b) ; the apocenter is reached at t ~ 10. 




X 



Fig. 13. — The plots refer to the cluster shown in 
Fig. 11 at t — 8.9. Upper panel: linear density 
as measured along rr-axis. Middle panel: aver- 
aged stellar velocity component (<v x >) along x. 
Lower: spatial derivative of <v x > (the horizontal 
line indicates the zero value). The position of the 
GC and of the 'ripple' are marked, respectively, by 
the dashed and dotted lines (galactic center is at 
x = 0). 
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Fig. 14. — As in Fig. 13 but referring to the time 
t = 10.5. The position of the GC is marked by 
the dashed line, while the 'clump' is at the dotted 
one. 
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